function result = DynesFit_nodal_gauss_fraction_Base (Param, meV) %Param(1) = delta(meV), Param(2)=gamma(meV), Param(3)=T(mK), Param(4) = slopeBG, Param(5) = offsetBG

delta0 = Param(1);
sigma = Param(5);
w = 1; % fraction of SC region
Temperature = 116.4; %Temperature of the sample (mK)
beta = 1/(8.61733E-5 .* Temperature);  % 8.61733E-5 is kb/e
xrange = -10/beta : 0.41/beta: +10/beta;  % Energy integration section for Fermi-Dirac distribution
yrange = 0 : 0.12 : 2*pi; %Gap integration across theta

if delta0-4*sigma < 0
zrange = 0 : 4*sigma/50 : delta0+4*sigma; % Delta integration section for gaussian distribution at E>0
zrange2 = delta0-4*sigma: 4*sigma/50 : 0; % Delta integration section for gaussian distribution at E<0

[XX,YY,delta] = meshgrid(xrange,yrange,zrange);
[XX2,delta2] = meshgrid(xrange,zrange2);

Gauss = 1./(sqrt(2*pi).*sigma).*exp(-(delta-delta0).^2/(2*sigma.^2)); %.* (1+erf(Param(7)*(delta-delta0)/sqrt(2))); %Gaussian distribution for E>0
Gauss2 = 1./(sqrt(2*pi).*sigma).*exp(-(delta2-delta0).^2/(2*sigma.^2)); %.* (1+erf(Param(7)*(delta2-delta0)/sqrt(2))); %Gaussian distribution for E<0

innerFtn1 = (1-w) .* beta .* exp(beta.*(xrange))./ ((exp(beta.*(xrange))+1.0).^2); % Normal state contribution
innerFtn2 = w .* Gauss2.* beta .* exp(beta.*(XX2))./ ((exp(beta.*(XX2))+1.0).^2); % State not SC yet contribution (E<0 in Gaussian)
innerFtn3 = @(v) w .* abs(real( (XX + v -1i.*Param(2)) ./ sqrt((XX + v -1i.*Param(2)).^2-(delta.*cos(YY)).^2))) .*Gauss.* beta .* exp(beta.* (XX))./ (2.*pi .*(exp(beta.* (XX))+1.0).^2);

result = arrayfun(@(v) ((Param(3)*v + Param(4)).* ( trapz(xrange,innerFtn1) + trapz(xrange,trapz(zrange2,innerFtn2,1)) + trapz(yrange,trapz(xrange,trapz(zrange,innerFtn3(v),3),2)))),meV);

elseif delta0-4*sigma > 0

zrange = delta0-4*sigma : 4*sigma/50 : delta0+4*sigma; % Delta integration section for gaussian distribution at E>0
[XX,YY,delta] = meshgrid(xrange,yrange,zrange);
Gauss = 1./(sqrt(2*pi).*sigma).*exp(-(delta-delta0).^2/(2*sigma.^2)); %.* (1+erf(Param(7)*(delta-delta0)/sqrt(2))); %Gaussian distribution for E>0
innerFtn1 = (1-w) .* beta .* exp(beta.*(xrange))./ ((exp(beta.*(xrange))+1.0).^2); % Normal state contribution
% innerFtn2 = w .* Gauss2.* beta .* exp(beta.*(XX2))./ ((exp(beta.*(XX2))+1.0).^2); % State not SC yet contribution (E<0 in Gaussian)
innerFtn3 = @(v) w .* abs(real( (XX + v -1i.*Param(2)) ./ sqrt((XX + v -1i.*Param(2)).^2-(delta.*cos(YY)).^2))) .*Gauss.* beta .* exp(beta.* (XX))./ (2.*pi .*(exp(beta.* (XX))+1.0).^2);

result = arrayfun(@(v) ((Param(3)*v + Param(4)).* ( trapz(xrange,innerFtn1) + trapz(yrange,trapz(xrange,trapz(zrange,innerFtn3(v),3),2)))),meV);
end

    
end